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Abstract 

A mathematical formulation and computational techniques are presented to describe optimal control and design 
strategies for the suppression of turbulent motions in the melt and the minimization of temperature gradients in the 
crystal in Czochralski crystal growth processes. The methodologies developed can be used to test control mechanisms, 
design parameters, and optimization objectives to determine their effectiveness in improving the processes. They can 
also be used to effect such improvements by systematically determining optimal values of the design parameters. The 
controls or design parameters considered include applied magnetic fields, temperature gradients along the side wall of 
the crucible, and crucible and crystal rotation rates. The results show that applied magnetic fields can be very effective in 
reducing velocity perturbations in the melt, while side wall temperature gradients are less effective and crucible and 
crystal rotation rates are ineffective. The results also show that applied magnetic field and crucible and rotation rates are 
ineffective in reducing temperature gradients in the crystal or in the melt. © 2002 Elsevier Science B.V. All rights 
reserved. 
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1. Introduction 

During the past several decades, three different 
types of crystal growth techniques have been 
developed, namely, vapor, melt, and solution [1]. 
Nevertheless, the melt growth technique, i.e., 
the Czochralski (Cz) process , has dominated the 
production of single crystals for most of the 
materials used in the microelectronics industry. 
However, in spite of the popularity of the Cz 
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method, there remain many limitations. For 
example, without an externally applied magnetic 
held, the melt motion in many commercial 
processes that use the Cz method is turbulent 
and involves large scale motions. FLeat flux 
fluctuations from the melt to the crystal produce 
a cycle of crystallization and remelting at the 
interface. This cycle produces dislocations and 
other microscopic defects in the crystal. Experi¬ 
ments on magnetic Czochralski (MCz) growth 
have shown that the temperature fluctuations 
and erratic doping striations are suppressed by 
the application of a magnetic field [2]. In fact, the 
most important benefit of the magnetic held is that 
it controls impurities and inhomogeneities at 
microscopic levels by producing better conditions 
in the vicinity of the melt/crystal interface. Never¬ 
theless, the use of MCz processes for the industrial 
production of crystals is still very limited because 
of its cost, as well as the limited understanding of 
when and how the magnetic field stabilizes the 
melt. Since MCz is expected to become the future 
technology for materials that are difficult to grow, 
particularly at high and ultra-high pressures, as 
well as the industry norm for the growth of large- 
size crystals, our current focus is on introducing 
magnetic effects into an existing Cz model. 

Here, we present a systematic approach to the 
control or optimization of the MCz crystal growth 
process. Even though the fundamentals of indivi¬ 
dual transport processes within the Cz process 
may appear to be simple, their interactions are 
rather complex. Thus, the notion of controlling the 
process is a nontrivial matter. In particular, it is 
seldom obvious which design parameters are 
effective for use in achieving a specific objective. 
Thus, our goal is to develop a methodology that 
can be used to test, through computational 
experiments, the effectiveness of design objectives 
and parameters. We also want this methodology to 
be useful for the actual design of improved crystal 
growth processes. To this end, we show that some 
design parameters, e.g., the strength of an applied 
magnetic field, can be very effective in improving 
some aspects of the MCz crystal growth process. 

We first investigate the effects on the melt flow 
velocity that occur as a result of various actions 
that can be applied, using design parameters such 


as the strength of an applied magnetic field, the 
temperature along the boundaries, and the crystal 
and crucible rotation rates. We also investigate the 
effect of these design parameters on the tempera¬ 
ture gradient in the crystal during the growth 
process. 

2. Mathematical model for the MCz process 

2.1. Governing e({nations 

A sketch of the physical domain used for our 
calculations is given in Fig. 1. The basic equations 
governing the melt in a MCz crystal growth 
process are a coupled magnetohydrodynamics 
(MHD) system. This system is composed of the 
conservation equations for fluid momentum, mass, 
and energy. Maxwell’s equations, and Ohm’s law 
for a medium in motion. The conservation 
equations for the melt, in primitive variables, are 
given by 

pQ^ + u- Vuj -/(Au + V/r = F b + F m , (1) 
V • u = 0, (2) 

^(f + u .vr)= K Ar, (3) 

where p, u,p, and T are the fluid density, velocity 
vector, pressure, and temperature, respectively. 



Fig. 1. Physical domain for Cz Si growth. 
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The parameters p, c p , and k appearing in the above 
equations are the fluid viscosity, specific heat 
capacity at constant pressure, and thermal con¬ 
ductivity, respectively. Fb is the buoyancy force, 
which is modeled in the Boussinesq approximation 
by [3,4] 

F b = -(p - p 0 )ge z = p Q P(T - Tf)ge z , (4) 

where e r denotes the unit vector in the z direction 
so that the gravitational acceleration is represented 
by the constant vector —ge z . In Eqs. (l)-(4), we are 
assuming that the variations in the density in the 
melt are negligible, i.e., the fluid density p = p Q is 
constant, except in the body force term Fb, where 
the density is taken to be p = p 0 (l — fi(T — Tf )). 
Here, p Q is the melt density, Tf is a reference 
temperature which we choose to be the melt 
freezing temperature, and /? is the thermal expan¬ 
sion coefficient. F m denotes the Lorentz force 
induced by the magnetic field on the moving fluid; 
it can dampen convection and can thus produces 
better flow conditions in the melt. 

For magnetohydrodynamic applications, Max¬ 
well’s equations and Ohm’s law are given by [5] 

j = <7(E + uxB), (5) 

VxB = /i e j, (6) 

V • B = 0, (7) 


VxE=--, (8) 

where a is electric conductivity of the melt, B is the 
magnetic induction, j is the induced current 
density, E is the electric field, and p e is the 
magnetic permeability of the melt. System (5) (8) 
may be combined to yield 

0B 1 

— =-AB + Vx(uxB) and V • B = 0. (9) 

at ap e 

We nondimensionalize utilizing the following 
scalings: length b, velocity v/b, pressure pv 2 /b 2 , 
time b 2 /v, and magnetic induction p e Ho, where 
v = p/p Q is the kinematic viscosity and Hq is the 
magnitude of the applied magnetic field. We set 
& = (T — Tf)/(T C — Tf), where T c is the constant 
temperature at the bottom crucible wall. Using 
these scalings, the nondimensionalized governing 


equations for flow and heat transfer in the melt are 
given by 

^ + (u • V)u + Vp - Au = Gr <9e z + F m , (10) 


V • u = 0, 


(ID 


^ + u . V 0 = 1 a0, 

at Pr 


( 12 ) 


where the components of F m in the x,r, and 9- 
directions are, respectively, given by 

F x = 0, F r = - ^-7-v, F e = -Ha 2 ^. (13) 

<7 m OX 


We solve the electric current stream function 
equation (ECSFE) instead of for the electric 
potential. The primary advantage of using the 
stream function equation is the simplicity in 
implementating boundary conditions. The bound¬ 
ary conditions for the potential equation are much 
more complicated on curvilinear walls. Then, 
(ECSFE) can be written as 


0 /er m 1 dW\ 0 /er m 1 0f'\ oH’ 
0.v\ a r 0x ) + 0x \ er r dr ) dx' 


The magnetic field is governed by 

0D j 

— = -—AB + Vx(uxB), V • B = 0 (15) 

dt Re m 

and the temperature in the crystal is governed by 

^7 = k s A0. (16) 

dt 

In these equations, Gr = g^b 3 (T c — T[)/v 2 , Pr = 
vcp/ic, Re m = p e cjv, Ha = p t Hgb\J a /p, and k s de¬ 
note the dimensionless Grashof number, Prandtl 
number, magnetic Reynolds number, Hartmann 
number, and nondimensional thermal diffusivity, 
respectively. 

Initial conditions can be specified on the melt 
velocity and temperature and the crystal 
temperature. 


2.2. Interface and boundary conditions 

We restrict attention to axially symmetric melt 
flows described in terms of cylindrical coordinates. 
Note that for computational reasons, an enclosing 
top surface has also been employed. 
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We implemented the boundary conditions easily 
by using the stream function equation. The 
boundary condition for Eq. (14) is 

X7'F = 0. (17) 


The evolving crystal/melt and free surface/melt 
interfaces have a profound effect on the quality of 
the grown crystal. Thus, in order to accurately 
capture this evolution, we model the melt region 
by 

Q m = {()',z)|0 

-o 0'Kz5$Z(r,0}, (18) 


where Z(r, t) is defined by 


Z(r,t) 


Z i()', t), 0 ^r^r*(t), 

Z 2 (r,t), r*(t)<r^b. 


(19) 


zo(r) is the dimensionless shape function for the 
bottom of the crucible. Z(r, /) is the dimensionless 
height of the moving interfaces, with Z\(r, t) 
representing the crystal/melt interface and Z 2 (r, t) 
the encapsulant/melt interface; r*(t) is the triple 
point where the melt, crystal, and encapsulant 
meet. The parameter b is thus chosen to be the 
distance from the axis of symmetry to the vertical 
wall of the crucible. Note that, in this setting, we 
are modeling the solidification of a pure substance 
with a fixed fusion temperature 7} and thus we 
assume that the solid and liquid phases are 
separated by a sharp interface given by s(z, r, t) = 
Zi(r, t) — z = 0. The energy balance at the crystal/ 
melt interface defines its movement; it is governed 
by [6] 

ezi 

dt 

Ste/ X's 00 s 00 n -/ 

Pr \7C m 072 077 


Wp(?) 



( 20 ) 


where u p (t) is the pull rate, Ste =C ps (2"' c - T[)/h{ is 
the Stefan number of the melt, /; p is the latent heat, 
and the subscripts s and m refer to crystal and melt 
regions. This equation is based on the assumption 
that the crystal and melt are not separated at the 
triple point. 

Denoting the free-surface position as z = 
Z 2 (r, t), the height of the free surface can be 


determined by solving 

0 2 Z 2 /0r 2 dZ 2 /dr 

(1 + (0Z 2 /0t-) 2 ) 3/2 + /*( 1 + (0Z 2 /07-) 2 ) 1/2 

= Bo(Z 2 — 2), (21) 

where Bo = pgb 2 /a & is the Bond number, c s is the 
surface tension, and the parameter A can be 
calculated from the melt conservation constraint 

f R ' Z ir dr + f Z 2 r dr = ^~, (22) 

Jo J Rr 2jt 

where V m (t) is the volume of the melt and 
Rr = r s /b is the radius ratio where r s is the radial 
co-ordinate of the crystal. 

There are two boundary conditions needed to 
solve for the shape of the encapsulant/melt inter¬ 
face. At the triple point, the meniscus is considered 
to be pinned to the edge of the crystal, and at the 
junction between the encapsulant and crucible 
wall, a 90-deg contact angle is assumed in view of 
the weak influence of the shape of the crucible 
meniscus on heat transfer. 

Additional boundary conditions for the pro¬ 
blem under consideration include the following: at 
the crystal/melt interface, 

u = Mint, v = Vint, w = Re s r, 0 = 0, (23) 

at the encapsulant/melt interface, 

0Z 2 0Z 2 8v du Ma 00 

dt V dr ’ 0/7 0t Pr 07? ’ 

-s— — Blfs(0fs ~ ©a), (24) 

7C S 07? 

at the bottom and side wall of the crucible, 
u = v = 0, w = Re c r, 0 = 0 W , (25) 

at the top and side wall of the crystal, 

-^=Bi s (0 s -0a), (26) 

07? 

where Re s = Q s b 2 /v and Re c = Q c b 2 /v are the 
crystal and crucible rotation Reynolds numbers, 
respectively, Q s and Q c are the rotation rates for 
the crystal and crucible, respectively, the sub¬ 
script fs refers to the free surface, Ma = 
( da s /dT)(T c — Tf)b/poi is the Marangoni number, 
0 W , 0 S , and 0 a are the dimensionless crucible wall 
temperature, crystal boundary temperature, and 
ambient temperatures, respectively, and 7?j nl and 
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Hint are the interface velocity components in the 2 
and r directions, respectively. Since the objective of 
our paper is to develop a control algorithm for 
quasi-state cases, the interface velocities will not be 
presented here. However, they are important in 
transient calculations; we calculated them but did 
not save them in simulations. 

Since we focused on controlling the melt flow 
with a magnetic field, we used the simplified 
radiation model which was developed by Prasad 
et al. [7]. It is based on the method of discrete 
exchange factor (DEF) [8], The radiative heat loss 
has been considered from the side wall of the 
crystal and the melt free surface with radiation 
Biot number Bi = eg {T 2 + T^)(T + T a )b/K s , 
where s is the emissivity and T a is the ambient 
temperature. 

3. Modeling MCz control 

There are many opportunities for introducing 
optimal control strategies into crystal growth 
processes. Ultimately, we want to apply control 
strategies in order to improve, either directly or 
indirectly, the quality of the crystal. Applying 
control strategies that directly affect crystal 
properties is very difficult. Therefore, as a first 
step in applying systematic optimal control strate¬ 
gies to crystal growth processes, we will apply 
control strategies that only indirectly affect crystal 
properties. 

Specifically, we will examine three different 
objectives for control along with, at first, two 
different control mechanisms. The first objective is 
to minimize the vorticity in the melt region; this is 
an effort directed at lessening large-scale turbulent 
motions inside the melt and to indirectly minimize 
oscillations in the melt-crystal interface. The 
second objective is to minimize the gradient of 
the temperature in the melt; this is an effort 
directed at lessening temperature variations in the 
melt. The third objective is to minimize the 
gradient of the temperature in the crystal; this is 
an effort directed at indirectly lessening the 
presence of residual stresses in the crystal. The 
two controls we will employ are the strength of the 
applied magnetic field and the temperature on the 


side wall of the crucible. For our nondimensiona- 
lized model, the strength of the applied field enters 
only through the Hartmann number; thus, the 
Hartmann number is one of our control para¬ 
meters. The temperature on the side wall of the 
crucible is assumed to be linear in the height 2 ; of 
course, the temperature at the bottom of the side 
wall should be the same as the temperature at the 
bottom of the crucible which is assumed to be 
constant, i.e., T = T c along the bottom of the 
crucible. Thus, on the side wall of the crucible we 
have that 

0(z)=\+t( Z ~ Zb ), (27) 

where 2 b = zo{b) and zl are the dimensionless 
heights at the bottom and top of the side wall of 
the crucible. Thus, the slope factor £ is a second 
control parameter at our disposal. Note that if 
£< — 1, then the temperature on the side wall at 
the top of the crucible will be lower than the 
melting temperature of the melt. This is, of course, 
unacceptable, since it would cause the melt to 
solidify at the crucible wall. Therefore, we impose 
the constraint £^ — 1 in our optimization studies. 
See Fig. 2 for the results of a simulation with 
£=—1.1; the figure gives level curves of the 
temperature in the melt and crystal. Note, for this 
value of £, the small region of solidified melt at the 
top of the side wall of the crucible. 


T 



Fig. 2. Nondimensionalized temperature in the melt and crystal 
for a linear crucible side wall temperature profile with { = —1.1. 
Note the small region of solidified melt adjacent to the top of 
the crucible wall. 
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We now give a mathematical description of the 
optimal control problems we consider. Here, (0, T) 
denotes the (nondimensionalized) time interval 
over which one wants to exercise control, f2 m 
denotes the melt region, and Q c denotes the crystal 
region. The first optimization problem is given by: 
final optimal values of Ha and £, such that 

[ [ |Vxu| 2 dQ dt (28) 

JO J 

is minimized, where a(t,r,z; Ha, £) satisfies govern¬ 
ing system (10)—(16) along with initial conditions 
and the boundary and interface conditions of 
Section 2.2. 

The second optimization problem is given by: 
fined optimal values Ha and £ such that 

.y 2 (Ha, j) = [ [ |VT| 2 dQ dt (29) 

J 0 J Qm 

is minimized, where T(t,r,z ; Ha, f) satisfies govern¬ 
ing system (10)—(16) along with initial conditions 
and the boundary and interface conditions of 

Section 2.2. 

The third optimization problem is given by: fined 
optimal values of Ha and c, such that 

Ha,£) = [ j \S7T\ 2 dQdt (30) 

Jo Jq c 

is minimized, where T(t,r,z ; Ha, f) satisfies govern¬ 
ing system (10)—(16) along with initial conditions 
and the boundary and interface conditions of 

Section 2.2. 

These represent nonconvex optimization 
problems. 

The solution algorithm for fluid flow calcula¬ 
tions in a generalized curvilinear coordinate 
system is basically similar to the SIMPLER 
algorithm [9] which consists of solving a pressure 
equation to obtain the pressure field and solving a 
pressure-correction equation to correct predicted 
velocities. We took initial velocity and temperature 
to be zero. 

We now lay the ground work for the computa¬ 
tional study of this optimization problem We 
apply the discrete Armijo gradient algorithm [10]. It 
should be noted that the nonconvexity of the 
functionals implies that the convergence of such an 
algorithm is conditional, i.e., it depends on having 


a good “initial guess”. We describe the algorithm 
in our current context; the search direction h, is a 
finite difference approximation to —V,/, with the 
parameter e controlling the precision of this 
approximation. 

The discrete Armijo gradient algorithm is a 
variable step size gradient method. It is particu¬ 
larly easy to implement since the step size selection 
does not require a one-dimensional line search. 
Convergence proofs under mild conditions on the 
functional are available [11]. In practice, we have 
found it to perform well, needing, on the average, 
approximately 15 iterations for satisfactory con¬ 
vergence but each iteration required five functional 
calculations in average. Of course, there are more 
sophisticated optimization methods available, e.g. 
trust-region quasi Newton methods. However, for 
our purposes, the Armijo algorithm was entirely 
adequate since our main goal in this study was to 
demonstrate the effectiveness of a systematic 
optimal control approach for improving processes 
such as Cz. 

Discrete Armijo Gradient Algorithm 


1. Choose a,j8e(0,1), ye(0, oo) and k*, k 0 e^f. 
Choose initial guesses for the control para¬ 
meters Ha and j. Set i = 0 and s = f k °. 

2. Determine the search direction h, having 
components 

, ,/(Ha + e, £) — ./(Ha, £) 

hi 1 =-(it) 

E 

and 


ha 


^(Ha,g + e)-^(Ha,Q 

e 


(32) 


3. Compute 

A _ ./(Ha + shn, £ + sha) — ./(Ha, £) 

23 i — 

£ 


(33) 


4. If Aj > 0, replace e by ySe, and go to Step 1. 
Else, use the subprocedure given below, which 
requires k*, to compute the step size f = [i k ‘, 
where k t e such that 

J(Ua + [l k di n ,c + fi k ‘h l2 ) 

-/(Ha,0<jd l H (34) 
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and 

/(Ha + p k ‘~ l ha,£, + p k ‘~ l ha) 

— ,/(Ha, £) > p kl ~ l aAj. (35) 

5. If Ha + khj ], + Aha) - ./(Ha, £) > —ye, re¬ 
place e by /?£, and go to Step 1. 

Else, replace Ha by Ha+AA/i, £ and c + hha, 
and i by i + 1; then, go to Step 1. 

Step size subprocedure 

1. If /' = 0, set k' = k*. Else, set k! = kj \. 

2. If kj = k! satisfied Eqs. (34) and (35), stop. 

3. If kj = k! satisfies Eq. (34) but not Eq. (35), 
replace k' by k' — 1 and go to Step 2. 

4. If kj = k' satisfied Eq. (35) but not Eq. (34), 
replace k! by k' + 1 and go to Step 2. 

3.1. Simulation code 

Over the past 25 years, transport phenomena in 
Cz crystal growth processes have been numerically 
simulated by many investigators; see, e.g.. Refs. 
[12-16], The simulation code used in our optimiza¬ 
tion studies is based on a high-resolution computer 
model developed in Refs. [6,17,18]. This model 
employs a multizone adaptive grid generation 
scheme and a curvilinear finite volume formula¬ 
tion. Once a grid has been generated for the flow 
domain, the control volume formulation begins by 
decomposing the domain into nonoverlapping 
control volumes such that there is one control 
volume surrounding each grid point. This compu¬ 
tational domain also allows for the coexistence of 
various materials in different phases with signifi¬ 
cantly different thermophysical and transport 
properties. Next, the governing MHD equations 
are integrated over each control volume. Piecewise 
profiles expressing the variation of state variables 
between the grid points are then used to evaluate 
the required integrals. The end result is the discrete 
equations containing the values of the state 
variables for a group of grid points. The appealing 
feature of the control volume formulation is that 
the integral conservation of mass, momentum, and 
energy are satisfied exactly over any group of 


control volumes and, of course, over the whole 
computational domain. For more details on the 
finite volume formulation, see, e.g.. Ref. [9], 


4. Numerical results 

In every case we studied, the optimal value of 
the side-wall temperature parameter turned out to 
be its limiting value £ = — 1. (Recall that we 
cannot allow £< — 1 since this would cause the 
melt to solidify at the side wall.) This is not a 
surprising result in optimization, i.e., the optimal 
value of a parameter is at the boundary of its 
admissibility set. As a result, we will not provide 
results for optimization studies with respect to 
Instead, we give results for £ = 0, i.e., a constant 
temperature profile, which is what was chosen in 
most past simulations, see, e.g., Refs. [6,7,19], and 
for £, = — 1, i.e., a linear temperature profile having 
as much slope as allowable, which provides the 
optimal results. 

The CPU time for solving the partial differential 
equations in 100 time steps with 102 x 50 mesh is 
about 23 min on Linux computer which had a 
600 MHz processor, Intel Pentium III with 
256 MBytes memory. 

For all the numerical experiments we describe 
here, we used 102 x 50 spacial grid and 100 time 
steps with a nondimensionalized time step 5t = 
0.0005; this corresponds to approximately 17 s. In 
fact, we did not see a significant change even we 
use 1000 time steps since we assume a quasi-state 
case. Otherwise, the melt drop will affect the heat 
transfer. For a silicon melt, we have the physical 
parameters a = 1.2xl0 6 S/m (electrical conductiv¬ 
ity), v = 3.0xl0~ 7 m 2 /s (kinematic viscosity), p = 
2420 kg/m 3 (density), and b = 0.1 m; see Table 1 
for more information. We assume that the 
magnetic field is in the axial direction, which is 
the most practical choice. The ambient tempera¬ 
ture is such that 0 = —6.00372; of course, 0 = 0 
and 1 corresponds to the freezing temperature of 
the melt and the constant temperature at the 
bottom of the crucible, respectively. 

The crystal radius has been specified in the 
simulations. The pulling rates can be determined 
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Table 1 

Thermophysical properties and parameters for Si melt and solid 


Properties of Si melt 


Density 

kg/m 3 

2420 

Thermal conductivity 

W/mK 

64 

Kinematic viscosity 

m 2 /s 

3.0 x 10~ 7 

Freezing temperature 

°C 

1410 

Surface tension 

N/m 

0.72 

Latent heat of solidification 

J/kg 

1.8 

Coefficient of thermal expansion 

/K 

1.41 x 10~ 5 

Specific heat capacity 

J/kgK 

1000 

Emissivity 


0.15 

Stefan number 


8.33 x 10~ 2 

Properties of Si solid 

Density 

kg/m 3 

2300 

Thermal conductivity 

W/mK 

22 

Specific heat capacity 

J/kgK 

1000 

Electrical conductivity 

S/m 

5.8 x 10 4 

Emissivity 


0.75 

Grashof number 


1.0 x 10 6 

Prandtl number 


15 x 10~ 3 

Stefan number 


4.17 x 10~ n 


Table 2 

Objective functional values vs. Hartmann number for different 
iterations of the discrete Armijo gradient algorithm (above the 
line) and for some additional values of the Hartmann number 
above the optimal value (bold face) determined by that 
algorithm (below the line). The imposed temperature along 
the side wall of the crucible is a constant ({ = 0) 


Iteration no. 

Hartmann no. 

Jf t 

1 

100.00000 

2121.056806 

2 

125.00000 

1286.909470 

3 

133.00000 

1143.171567 

4 

134.66505 

1117.664065 

5 

134.98452 

1113.226481 

6 

135.09339 

1112.402257 

7 

135.11556 

1112.379587 

8 

135.11132 

1112.378560 

9 

135.11178 

1112.378542 

— 

135.11200 

1112.378545 

— 

136.12222 

1169.954222 

— 

138.00000 

1326.942849 

— 

140.00000 

1532.767010 

— 

150.00000 

4043.390458 


based on the given radius. The pulling rates are 
small compared to the melt velocity. 

4.1. Minimization of the vorticity functional J\ 

4.1.1. Constant temperature wall profile 

For the choice £, = 0, a constant side-wall 
temperature, we initially chose the parameters for 
the discrete Armijo gradient algorithm as follows: 
e = 0.01, k = 6, [l k = e, i.e., 0.46, a = 0.5, y = 5, 

and FIa= 100. Table 2 gives, for each optimization 
iteration, the value of the Hartmann number and 
the corresponding value of the functional ,/j. The 
optimization algorithm converged after nine itera¬ 
tions; the optimal value of the Hartmann number 
for this case is Ha = 135.1117816 which, for our 
choice of physical parameters, corresponds to an 
applied magnetic field of 0.0331 T; u= 113.3, 
velocity in axial direction, corresponds to 
0.0003m/s; Re cru =—500, crucible rotation rate, 
corresponds to —0.0151/s. We also computed the 
value of the functional for values of the Hartmann 
number above the optimal value produced by the 
optimizer in order to show that the optimizer did 


indeed find a minimum of the functional. These 
results are also given in Table 2. In Fig. 3, we 
provide plots of the values of the logarithm of the 
functional and of the Hartmann number vs. the 
iteration number; these clearly show the 
convergence of the iterative method. In that 
figure, we also provide a plot of the value of 
the logarithm of the functional vs. the Hartmann 
number, where the latter is chosen to be the 
values visited during the optimization iteration 
and some additional values chosen after the 
iteration is completed. We clearly see that the 
optimizer did indeed find a minimum of the 
functional. In Figs. 4 and 5, we provide plots of 
the velocity vector and of its magnitude for no 
control applied, i.e., no applied magnetic filed 
(Ha = 0), and for the optimal magnetic field 
(Ha = 135.11178). Notice that convection in the 
melt is significantly reduced for the optimal value 
of Ha. In fact, the maximum magnitude of the 
velocity decreases from 694.547 for Ha = 0 to 
27.6721 for the optimal value of Ha. In Figs. 4 and 
5, we also provide contour plots of the tempera¬ 
ture in the melt and the crystal. 
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Fig. 3. Functional and Hartmann number values for a constant crucible side wall temperature profile (£ — 0). (a) Value of the 
functional vs. optimization iter. Number as determined by the optimization process (first nine points) and for some additional 
values of the Ha number (last five points), (b) Ha number vs. optimization iter. Number as determined by the optimization process 
(first nine points) and some additional values chosen to determine if the optimizer truly found a minimum (last five points), (c) Value of 
the logarithm of the functional vs. Hartmann number showing the existence of a minimum. The values of the Ha number for the 
points plotted to the left of the minimum were determined by the optimization process; those to the right were chosen so as to show the 
minimum. 
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4.1.2. Optimal temperature wall profile 

For the choice £ = — 1, the optimal linear side- 
wall temperature, we initially chose the parameters 
for the discrete Armijo gradient algorithm as 
follows: s = 0.001, k = 6, fl k = e, i.e., fix 0.32, a = 
0.5, y = 7, and Ha= 10. After five iterations there 
was a need to change a to 0.86 and /? to 0.35. 
Table 3 and Figs. 6-8 provide, for £ = —1, the 
corresponding information that was provided for 
£ = 0 in Table 2 and Figs. 3-5, respectively. The 
optimization algorithm converged after 21 itera¬ 
tions; the optimal value of the Hartmann number 
was now found to be Ha = 134.66505. Again, the 


convection in the melt is significantly reduced for 
the optimal value of Ha. In the Ha = 0 case, a large 
convective flow was observed in the melt and more 
specifically near the interface; for Ha= 134.66505, 
which was obtained by optimization, we see a 
greatly reduced convective flow. For t = —1, the 
maximum magnitude of the velocity decreases 
from 507.8 for Ha = 0 to 10.55 for the 
optimal value of Ha in Table 4. Of particularly 
interest are the results along the crystal interface; 
there the norm of velocity was reduced from 66 for 
Ha = 0 to 10.6 for the optimal value 
Ha = 134.66505. 
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Fig. 4. Temperatures in the melt and crystal and velocities in the melt for a constant crucible side wall temperature profile ({ = 0) and 
Fla = 0 (no magnetic control applied); (a) the left half gives the contours of the speed and the right half gives the velocity vector, (b) the 
left half gives the temperature and the right half gives the velocity vector. 



Fig. 5. Temperatures in the melt and crystal and velocities in the melt for a constant crucible side wall temperature profile ({ = 0) and 
Fla= 135.11178 (the optimal value determined by the optimizer for the functional ,f\). (a) The left half gives Ihe contours of the speed 
and the right half gives the velocity vector, (b) the left half gives the temperature and the right half gives the velocity vector. 


4.1.3. Comparison of constant and linear wall 
temperature profiles 

From Table 4 and Figs. 9 and 10, it is clear that 
the linear temperature profile (27) with £ = — 1 
along the crucible wall yields better results, e.g., a 
greater suppression of convective perturbations in 
the melt. For example, from Table 4, we see that 
for Ha = 0 there is a reduction of the maximum 
speed of flow perturbations from 694.5 for the 
constant profile to 507.8 for the linear profile. For 


the corresponding optimal Hartmann numbers, 
the reduction is from 27.67 to 10.55. Plots of the 
velocity, stream function, and temperature for the 
constant and linear temperature crucible wall 
temperature profiles are given in Figs. 9 and 10. 
It is clear from the stream function contour plots 
that the stream function gradient has been 
reduced, which shows that convection in the melt 
has also been reduced; see Figs. 9a and 10a. From 
Table 4, we see that, for the corresponding optimal 
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Table 3 

Objective functional values vs. Hartmann number for different 
iterations of the discrete Armijo gradient algorithm (above the 
line) and for some additional values of the Hartmann number 
above the optimal value (bold face) determined by that 
algorithm (below the line). The imposed linear temperature 
profile along the side wall of the crucible has maximum 
allowable slope (£ = — 1) 


Iteration no. 

Hartmann no. 


1 

10.00000 

49455.851095 

2 

34.30591 

7510.535628 

3 

40.47653 

4663.203345 

4 

51.09596 

2378.748929 

5 

64.24992 

1279.573829 

6 

80.41390 

747.534616 

7 

100.96154 

476.158066 

8 

109.16027 

418.440560 

9 

115.15040 

386.126853 

10 

119.97603 

364.650641 

11 

124.05711 

349.064600 

12 

127.61256 

337.116627 

13 

130.77333 

327.605411 

14 

133.62488 

319.821284 

15 

134.53541 

317.518290 

16 

134.62291 

317.367216 

17 

134.63926 

317.349610 

18 

134.65084 

317.340238 

19 

134.65845 

317.335694 

20 

134.66313 

317.333594 

21 

134.66505 

317.332895 

_ 

136.22639 

418.186535 

— 

138.92060 

691.996673 


values of the Hartmann number, the spread in the 
stream function is reduced from 1.86 for the 
constant temperature profile to 0.66 for the linear 
temperature profile. 

4.1.4. Effects of variations in the Hartmann number 
So far we have examined the effect of two design 
parameters, the Hartmann number corresponding 
to an applied axial magnetic field and the slope in 
the linear temperature along the crucible wall. We 
have seen that optimization with respect to either 
of these can result in a reduction in the perturba¬ 
tion velocity and the vorticity in the melt. We have 
also seen that variations in the Hartmann number 
can be much more effect than variations in the 
slope of the wall temperature profile. For example. 


at Ha = 0, going from constant wall temperature 
profile to the linear profile with maximum allow¬ 
able slope effected a reduction in the velocity norm 
from 694.5 to 507.8. However, for each of these 
extremes in the slope of the wall temperature 
profile, much greater reductions can be effected by 
optimizing with respect to the Hartmann number. 
For the constant temperature profile (£ = 0), the 
perturbation velocity norm was reduced from 
694.5 at Ha = 0 to 27.67 at the optimal value 
Ha= 135.11178 and for the best linear temperature 
profile (c = —1), the perturbation velocity norm 
was reduced from 507.8 at Ha = 0 to 10.55 at the 
optimal value Ha= 134.66505. 

The effects of variations in the Hartmann 
number for the optimal linear wall temperature 
profile is illustrated in the bottom plots of Figs. 7 
and 8 where the temperature is plotted for Ha = 0 
and the optimal value Ha = 134.66505. From these 
plots, it is clear that the magnetic field has a 
smoothing and stabilizing effect on the tempera¬ 
ture distribution in the melt. 

4.1.5. Effects of other design parameters 

There are a number of other design parameters 
that may be used in attempting to reduce velocity 
and vorticity perturbations in the melt. For 
example, two that have been suggested are the 
crucible and crystal rotation rates. Before using 
these parameters in our sophisticated optimization 
methodology, we examined how variations in their 
values affected the velocity perturbations in the 
melt. Variations in the Hartmann number have a 
large effect on the maximum and minimum values 
of the velocity components. However, variations in 
the crucible and crystal rotation rates seem to have 
no effect on the velocity perturbations. Thus, these 
are not effective design parameters to use if one 
wants to suppress velocity perturbations in the 
melt. 

4.2. Minimization of temperature gradient 
functionals J 2 in melt and in crystal near 
interface 

We have observed how some design parameters 
(e.g., the Hartmann number) are very effective in 
reducing velocity perturbations in the melt, others 
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Fig. 6. Functional and Flartmann number values for (lie optimal linear crucible side wall temperature profile ({ = —1). (a) Value of the 
functional vs. optimization iteration number as determined by the optimization process (first 21 points) and for some additional 
values of the Fla number (last two points), (b) Fla number vs. optimization iteration number as determined by the optimization process 
(first 20 points) and some additional values chosen to determine if Ihe optimizer truly found a minimum (last two points), (c) Value of 
Ihe functional J\ vs. Fla number showing the existence of a minimum. The values of the Ha number for the points plotted to Ihe left of 
Ihe minimum were determined by the optimization process; those to the right were chosen so as to show the minimum. 


are mildly effective (e.g., the slope of a linear 
temperature profile along the wall of the crucible), 
while others are very ineffective (e.g., the crucible 
and crystal rotation rates.) We now want to see 
how effective these parameters are in reducing 
temperature gradients in the melt, i.e., for the 
minimization of the functional <# 2 , and in the 
crystal, i.e., for the minimization of the functional 
J?,. Again, before using these parameters in our 
sophisticated optimization methodology, we ex¬ 
amined how variations in the values of the design 
parameters affect the functionals $2 and ,/ 3 . A 
sampling of results is given in Table 5. We see that 


variations in the Hartmann number and the 
crucible and crystal rotation rates have almost 
no appreciable effect on the values of J 2 and 
which are root-mean-square measures of the 
temperature gradient. However, variations in the 
control parameters can effect some reductions in 
the maximum value of the temperature gradient 
and in the values at specific points and subregions 
such as in the triple point where the crystal, melt 
and encapsulant meet and the areas near the 
crystal/melt interface. Therefore, optimization of 
functionals such as J 2 and ,/ 3 is ineffective for 
reducing temperature gradients in the crystal. This 
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Fig. 7. Temperatures in the melt and crystal and velocities in the melt for optimal linear crucible side wall temperature profile ({ = — 1) 
and Ha = 0 (no magnetic control applied) (a) the left half gives the contours of the speed and the right half gives the velocity vector, 
(b) the left half gives the temperature and the right half gives the velocity vector. 



Fig. 8. Temperature in the melt and crystal and velocities in the melt for optimal linear crucible side wall temperature profile (q = — 1) 
and Ha= 134.66505 (the optimal value determined by the optimizer for the functional ,f\). (a) The left half gives the contours of the 
speed and the right half gives the velocity vector, (b) the left half gives the temperature and the right half gives the velocity vector. 


Table 4 

Maximum and minimum melt flow temperature and stream function values and maximum speed for no applied magnetic field and for 
Ihe optimal values (with respect to the functional ,/i) of the applied field strength and for constant ({ = 0) and best linear crucible wall 
temperature profiles ({ = — 1) 


Flartmann no. 

Temperature 



Stream function 

Max speed 

Wall profile 

Min 

Max 

Min 

Max 


0 

Constant 

-6.003 

1.0 

-86.81 

0.59 

694.5 

0 

Linear 

-6.003 

1.484 

-64.24 

2.92 

507.8 

135.11178 

Constant 

-6.003 

1.0 

-1.864 

0.0 

27.67 

134.66505 

Linear 

-6.003 

1.484 

-0.665 

0.309 

10.55 
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Fig. 9. Stream function and velocity in the melt and temperature in the crystal and melt for the optimal values of the Hartmann 
number determined by the minimization of the functional and { = 0 (a constant crucible side wall temperature profile), (a) Ihe left 
half gives the contours of the stream function and the right half gives the velocity vector, (b) the left half gives the temperature and the 
right half the velocity vector. 



(a) (b) 


Fig. 10. Stream function and velocity in the mell and temperature in the crystal and melt for the optimal values of the Hartmann 
number determined by the minimization of the functional ,f\ and £ = — 1 (a constant crucible side wall temperature profile), (a) the left 
half gives the contours of the stream function and the right half gives the velocity vector, (b) the left half gives the temperature and the 
right half the velocity vector. 


is due to the small value of the Prandtl number. 
Recall that the Prandtl number is the ratio of the 
momentum and thermal diffusivities, i.e., it is a 
measure of the relative size of the momentum 
transport and conductive heat transfer effects. We 


should point out that in the crystal region there 
were other local maxima of the temperature 
gradient in regions for away from the melt region; 
these values were very insensitive to changes in 
parameter values. 
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Table 5 

Maximum temperature gradient in melt and crystal and functional values in melt and crystal for several values of the Hartman 
number and of the crystal and crucible rotation Reynolds numbers 


Hartmann no. Crucible rotation rate 

Crystal rotation rate I^Cmax 

In melt 


In crystal 

Functional 

In melt 

value 

In crystal 

0 0 

1000 

39.3172 


2.21426 

0.06796 

15.25665 

-500 

0 

38.5130 


2.2693 

0.06631 

15.25669 

-500 

2000 

36.1654 


2.2252 

0.06692 

15.25668 

-1000 

1000 

39.5406 


2.3522 

0.06645 

15.25668 

-3000 

1000 

38.7061 


2.3035 

0.06669 

15.25667 

30 0 

1000 

35.8517 


2.0597 

0.06479 

15.25671 

-500 

0 

37.9086 


2.2316 

0.06484 

15.25670 

-500 

2000 

36.1654 


2.2252 

0.06692 

15.25668 

-1000 

1000 

37.8650 


2.2155 

0.06478 

15.25671 

-3000 

1000 

37.1574 


2.1715 

0.06498 

15.25669 

134.66504 0 

1000 

35.3139 


2.0896 

0.06439 

15.25674 

-500 

0 

35.9296 


2.1215 

0.06439 

15.25673 

-500 

2000 

33.9489 


2.0307 

0.06440 

15.25674 

-1000 

1000 

35.4594 


2.0883 

0.06439 

15.25673 

-3000 

1000 

34.8467 


2.0454 

0.06441 

15.25672 

150 0 

1000 

35.2913 


2.0884 

0.06440 

15.25674 

-500 

0 

35.8170 


2.1156 

0.06439 

15.25673 

-500 

2000 

34.0652 


2.0344 

0.06440 

15.25674 

-1000 

1000 

35.3645 


2.0854 

0.06439 

15.25673 

-3000 

1000 

34.8196 


2.0473 

0.06441 

15.25672 

199 0 

1000 

35.2751 


2.0852 

0.06446 

15.25674 

-500 

0 

35.4982 


2.0974 

0.06446 

15.25673 

-500 

2000 

34.4264 


2.0474 

0.06446 

15.25674 

-1000 

1000 

35.2657 


2.0804 

0.06446 

15.25673 

-3000 

1000 

34.5574 


2.0375 

0.06446 

15.25672 

5. Concluding remarks 


have a drastic 

effect on 

the success of an 



optimization 

or 

control strategy. For example, 

In this paper, our goal has been to demonstrate 

we have shown that optimization with respect to 

some of the techniques being used to develop an 

the slope of a 

linear temperature profile along the 

integrated approach for intelligent 

modelling. 

crucible wall 

has a mild 

effect on the size of 


design and control of crystal growth processes. velocity and vorticity perturbations in the melt. On 

This is a challenging task; however, the potential the other hand, drastic reductions can be effected 

benefits are large. Our current thrust has been on by optimizing with respect to the strength of an 

developing control and optimization algorithms applied axial magnetic field and a negligible effect 

that can be implemented into existing codes for the is realized by optimizing with respect to crucible 

MCz crystal growth process. Carrying out this and crystal rotation rates. This illustrates how, for 

goal has involved an integration between the a specific design objective, optimization with 

experiments, modeling, simulations and control respect to some design parameters is much more 

experts and results. effective than with respect to others. 

In our studies we have shown how the choices of We have also shown that although optimization 


design parameters and objective functional can with respect to the Hartmann number can 
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drastically reduce the size of velocity and vorticity 
perturbations in the melt, it has negligible effect on 
the size of the temperature gradient in the crystal. 
This illustrates how the same design parameter has 
little influence on some functionals and great 
influence on others. 

Current and future work includes the refinement 
and further development of practical control and 
optimization strategies; the integration of our 
results and algorithms into the latest version of 
MASTRAPP3d which simulates three dimensional 
MCz crystal growth processes; and the experi¬ 
mental verification and system integration of the 
control and optimization results obtained through 
computations. 
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